% clear all;
% clc;

% ----------------------- Importing Data ----------------------------------

patnum=1;
ptag=['p',num2str(patnum),'_'];

noaxons=200;
epidur=0;
th_elec=0; % 0, -10, and -20
patgeom=load('patient_cordgeom.txt');

if(epidur==1)
    loctag='epi_';
    ftag=['d',strrep(num2str(patgeom(patnum,8)),'.','p'),'_'];
else
    loctag='sub_';
    ftag=['d',strrep(num2str(patgeom(patnum,9)),'.','p'),'_'];
end
thtag=['t',strrep(num2str(th_elec),'-','n')];
dcpre=[ptag,'neg_',loctag,'bp_dc_act_',ftag,thtag];
drpre=[ptag,'neg_',loctag,'bp_dr_act_',ftag,thtag];

cd(['patient',num2str(patnum),'_potentials']);

uni_cur=load([ptag,loctag,'bp_current_',ftag,thtag,'.txt']);

cd ..;

cd(['patient',num2str(patnum),'_thresholds']);

dc_th=abs( load([dcpre,'.txt']) );
dc_th(:,1)=dc_th(:,1)/(2/uni_cur)*1e3;
dcbyderm = zeros(10,10);
dcbyderm(:) = dc_th(1:(noaxons/2),1);
dr_th=abs( load([drpre,'.txt']) );
dr_th(:,1)=dr_th(:,1)/(2/uni_cur)*1e3;
cd ..;

aa = 1:10;
bb = (aa'*ones(1,10))';
xx = bb(:);
yy = dcbyderm(:);

numax=10;
for ii=1:numax
    
    ko=numax*(ii-1)+1;
    kf=numax*ii;
    
    for jj=1:(numax-1)
        
        kk=numax*(ii-1)+jj;
        yycomp=[zeros(jj,1);yy((kk+1):kf)];

        ck1=aa(yy(kk)==yycomp);
        
        if(any(ck1))
            kk_p=numax*(ii-1)+ck1; 
            xx(kk_p)=xx(kk_p)+(1:length(ck1))*(2/numax);
        end
        
    end
end

% rtt1 = 4/0.011493/1e3;
% ratp = 4/0.005924/1e3;
% yy = yy/rtt1;

% ---------------------- Plotting Results ---------------------------------

thtag = ['Lateral Deviation = ',num2str(th_elec),'^o'];
t8_vertlv  = {'T9','T10/11','T12','L1','L2/3','L3/4','L5','S1/2','S2/3','S4/5'};

figure;
hold on;
plot(xx,yy,'k.','MarkerSize',30);
plot([0,11],[min(dr_th(:,1)),min(dr_th(:,1))],'k--','LineWidth',4);
ttl = ['Patient ',num2str(patnum),', ',thtag];
title(ttl,'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',36);
hold off;
xlim([0,11]);
ylim([0.15,2]);
set(gca,'XTick',[],'YScale','Log','YTick',[0,0.1,0.2,0.5,1,1.5,2],'FontSize',36);
%set(gca,'YScale','Log');
for jj = 1:length(t8_vertlv)
    nolet = length(t8_vertlv{jj});
    hjj=text(jj-0.175*nolet/2,0.11,t8_vertlv{jj},'FontSize',30);
    set(hjj,'rotation',60);    
end
legend('DC fibers','DC_0','Location','NE');
axis square;
